Network formation of tissue cells via preferential attraction to elongated structures 
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Vascular and non-vascular cells often form an interconnected network in vitro, similar to the early 
vascular bed of warm blooded embryos. Our time-lapse recordings show that the network forms by 
extending sprouts, i.e., multicellular linear segments. To explain the emergence of such structures, 
we propose a simple model of preferential attraction to stretched cells. Numerical simulations 
reveal that the model evolves into a quasi-stationary pattern containing linear segments, which 
interconnect above the critical volume fraction of 0.2. In the quasi-stationary state the generation 
of new branches offset the coarsening driven by surface tension. In agreement with empirical data, 
the characteristic size of the resulting polygonal pattern is density-independent within a wide range 
of volume fractions. 



Embryogenesis, the shaping of tissues and organs, is 
a long standing challenge of science. Recent applica- 
tion of ideas originating from non-equilibrium statisti- 
cal physics has proved to be fruitful for understand- 
ing emergent properties such as pattern formation, or 
biomolecular network function during embryo develop- 
ment [1]. Formation of the primordial vascular bed of 
warm blooded vertebrates is arguably one of the best 
examples for an emergent phenomenon in embryonic de- 
velopment [2]. During vasculogenesis, well before the on- 
set of circulation, hundreds of essentially identical vas- 
cular endothelial cells create a polygonal network within 
a simple, sheet-like anatomical environment (Fig. [T^,b). 
The network is highly variable among individuals within 
certain statistical constraints relatively uniform in 
morphology, and unlike insect segmentation, there is no 
evidence for direct genetic control of vascular segment 
positions. 

The formation of linear segments and their intercon- 
nected network is not restricted to vascular endothelial 
cells [H]. In particular, we demonstrate in Fig. [I] (c) and 
(d) that non-vascular, glia- or muscle-related cells also 
exhibit linear structures when grown on a rigid plastic 
tissue culture substrate with continuously shaken culture 
medium. Depending on the cell density, the linear seg- 
ments can merge and form an irregular network. 

Recent in vivo observations of vasculogenesis indicated 
that early vascular network formation includes sprout- 
ing, the extension of linear segments containing multi- 
ple cells This process is markedly different from the 
gradual coarsening of an initially uniform density field, 
and its possible arrest, characteristic for colloid gels (see, 
e.g. 0). To determine the mechanism of sprout forma- 
tion, we recorded the temporal development of the struc- 
tures shown in Fig. [Hb-d) with computer controlled mi- 
croscopy [B[ . In vitro culture conditions yield sufficiently 
high resolution to trace the motion of individual cells 
during the patterning process. As shown in Fig. [TJe,f), 
sprout expansion involves cell motility guided by adja- 
cent projections of other cells or elongated multicellular 



structures. 

A number of theoretical models were formulated to un- 
derstand the self-organization of vascular networks. The 
chemo-mechanical mechanism assumed cells to exert me- 
chanical stress on the underlying substrate, and the re- 
sulting stress to guide cell motility [l^, 11 1. A recent 



model of Gamba, Serini et al assumed chemoattractant 



signaling [12.1131]. While the suggested chemoattractant, 
VEGFi65, is unlikely to behave in the predicted way (see 
the discussion in [14[), patterning guided by an unspec- 
ified chemoattractant continues to serve as the basis of 
biologically plausible models including cell adhesiveness 
and finite cell size [3]. 

Both the chemo-mechanical and chemoattractant 
mechanisms may be biologically relevant under certain 
circumstances. We argue, however, that neither can ex- 
plain sprout formation seen in Fig. [1] In both models 
the pattern were reported to arise in a gradual coarsening 
process in which segments are eliminated and small holes 
are replaced with larger ones. Moreover, neither mech- 
anism is expected to operate within our in vitro exper- 
imental setup: The rigid substrate excludes the chemo- 
mechanical mechanism. A specific chemotactic response 
is unlikely to be shared by such a variety of cell types. 
Finally, convection currents in the culture medium, gen- 
erated by temperature inhomogeneities within the incu- 
bator and the vibrations of microscope stage motion, are 
expected to hamper the maintenance of concentration 
gradients, or impose a strong directional bias upon the 
chemotaxis-related cell movements. 

To measure the convection currents close to the cul- 
ture surface, we immersed 0.5 /xm diameter latex beads 
(Sigma) into the medium. Bead motion was recorded 
within a 20 fim thick volume above the culture sur- 
face, delimited by the field depth of the 10X microscope 
objective. As a representative sample (Fig. Id, inset) 
demonstrates, in our experimental setup convection cur- 
rents were sustained for hours with speeds exceeding 100 
/im/h, an order of magnitude larger than the median cell 
speed. 
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FIG. 1: Tissue cell networks, (a) Vascular endothelial cells 
during vasculogenesis, visualized in a bird embryo according 
to [f| . (b) Mouse vascular endothelial cells within an in vitro 
explant [Q], after 24 h of culturing. (c) Astroglia-related rat 
C6 cells and (d) muscle-related mouse C2C12 cells cultured on 
a rigid tissue culture substrate as described in Q. The inset 
depicts the motion of latex beads at the tissue culture surface. 
To measure convection currents, 10 /im wide and 250 /im long 
rectangles were selected in the images, directed parallel to the 
flow. For each frame, pixel intensities were averaged across 
the rectangle, resulting in the horizontal lines of the inset, 
(e) C6 cells (arrows) migrate along extended projections of 
adjacent cells (arrowheads), (f) Endothelial cells (arrows) 
move along a vascular sprout, visualized in the mouse explant 
by differential interference contrast microscopy. 



On the basis of our empirical observations, in this let- 
ter we propose that cellular sprouting employs a quite 
generic mechanism, the preferential attraction to elon- 
gated structures. While the molecular basis of such a 
behavior is unknown, one may conjecture two biologi- 
cally feasible scenarios, (i) The tendency of cells to align 
with one another was analyzed in detail in Myxobacteria 
(see, e.g. [HI), and similar mechanisms may also operate 
in animal cells, (ii) Cells in elongated structures are as- 
sumed to be under mechanical tension, and strained cells 
can have a stiffer cytoskeleton [l6j]. Cells are able to re- 
spond to variations in extracellular matrix stiffness 



FIG. 2: Network formation in the model. Randomly placed 
N — 500 particles assemble into linear structures, detectable 
already within 30' (a). At sufficiently high particle density a 
characteristic pattern size develops in five hours (b) with a 
combination of sprouting (branch extension) and coarsenen- 
ing (merger of adjacent branches). Connected dots represent 
Voronoi neighbor particles. Darkening gray levels indicate in- 
creasing local anisotropy. The simulation covered an area of 
L — 0.7mm. 



is also feasible. To assess the collective behavior of cells 
exhibiting the proposed preferential attraction property, 
we studied the following simple model in which individual 
cells are represented as particles. 

Cell motility is often approximated as a persistent dif- 
fusion process mm, 

where the velocity Vk of cell k is 
described by the Langevin equation 



dvk 
~~dt 



Vk/r+VDZk+Mk 



(1) 



where r is the persistence time of motion, D is a dif- 
fusion parameter, and £ is an uncorrelated white noise: 
(0 = and (&(*)&(*')) = S H 6(t - t'). Term M is a 
deterministic bias, representing interaction with the en- 
vironment (specified below). While Eq.((T]) describes the 
motion of a Brownian particle at finite temperatures, 
animal cell motility is driven by complicated molecu- 
lar machinery, and it is not thermal fluctuation driven. 
Thus, parameters r and D depend substantially on cell 
type and molecular state. Measurements performed with 
non-interacting endothelial cells and fibroblasts resulted 
t and D values in the 01 - 5 h and 100-2000 ^m 2 /li 3 
range, respectively 
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Interactions among mobile agents are usually modelled 
as a sum of pair interactions [2l[ . In this spirit, we factor 
M into 



Mi, 



E 



Xk 



dk 



[h{d kj )+Wih(d kj )\, (2) 



and an analogous mechanotaxis utilizing cell-cell contacts 



where the sum is taken over the Voronoi neighbors of 
particle k, and dkj — \xj — x k \- The soft-core repul- 
sion fi ensures that model cells are impenetrable. The 
range of repulsion is the size Ri of the organelle-packed 
region around the cell nucleus. The preferential attrac- 
tion response is incorporated in the f% term. Cells are 
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expected to explore their surroundings with protrusions, 
and respond when protrusions contact elongated struc- 
tures. Filopodia typically extend from R 2 , the cell sur- 
face (Ri < R 2 ) 1 to a maximal distance of R. Thus, 
fi(d) = for d > Ri and f 2 (d) = for d < R 2 and d > R. 
Based on Fig.Q]we estimate R\ = lOyum, R 2 = SOfim and 
R = 40yum. These values, however, can vary by at least 
a factor of two, depending on the cell types and experi- 
mental conditions. 

For the sake of simplicity, cell shape is not explicitly 
resolved in the model. Therefore, elongation or local 
anisotropy is inferred from the configuration of particles. 
To represent an attraction to cells within anisotropic 
structures, the weights Wk are constructed as 



Wk 



1 



£ < 

{j:d jk <R} 



(3) 



where the sum is taken over all particles that are 
within a circle of radius R around particle k. The angle 
between x k — Xj and an arbitrary reference direction is 
denoted by tpjk- Thus, w — for particles in an isotropic 
environment and w — 1 for particles in a highly elon- 
gated, linear configuration. Non-uniform weights result 
in asymmetric pair-interactions, which is feasible as M 
represents a bias in cell activity instead of physical forces. 

Eqs (l)-(3) were studied by Euler integration with 
0.05h long time steps within a rectangular area of size L. 
Random initial and periodic boundary conditions were 
applied. Parameter values r = 0.5h and D = 100/im 2 /h 3 
were chosen to represent cell motility. There is little em- 
pirical guidance on the choice of functions fi and f 2 . 
We tested (i) a linear form with R x = R 2 = 15/zm as 
fi(d) + f 2 (d) ~ d — i?i and (ii) a Hertz- type repulsion 
fi(d) = — A(R\ — d) n with a zonal, distance-independent 
attraction f 2 (d) — B for R 2 < d < R, where the Hertz 
exponent n, A and B are parameters. Irrespective of the 
functional forms (i) and (ii), as well as for both n = 1 
and n = 1.5 a parameter regime exists in which sprout- 
like linear structures form during a biologically feasible 
time period. In Fig 2. and subsequently we demon- 
strate simulation results of model (ii) with parameters 
n = 1, A = 160 hr 2 and B = 130 ^m/h 2 . These val- 
ues represent a strong response to external cues: the ra- 
tio of the directed and random velocity components is 
Br I sj Dt f« 3. As a comparison, a similar measure for 
chemotactic response of endothelial cells was found larger 
than one [l9j ]. 

The time development of the system was followed by 
the ensemble-averaged binned power spectrum S(q) — 
(I J2j ex P( — 27Ti<fa?j)| 2 ), where each component of q is an 
integer multiple of 1/L. The (...) average is taken over 
configurations obtained with different noise realizations 
and over all scatter vectors q of length q. As Fig. 3a 
demonstrates, after an initial coarsening regime the pat- 
tern reaches a quasi-stationary state in which the gener- 
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FIG. 3: Power spectrum S(q) of the particle configurations, 
at various time points and model parameters. The dotted, 
dashed and solid curves indicated by the arrows were obtained 
at t =25, 75, 125 and 500h, respectively. The asymmetric 
model (a) reaches a steady state while the symmetric model 
(b) continues to coarsen into droplets (inset). The shifted 
curves in (a) demonstrate that q* is independent of the system 
size and cell density in the stationary phase. For low cell 
densities the peak falls off as 1/q, in good agreement with 
data obtained for Fig. (filled symbols). The solid lines 
represent power-law decays with exponents —1 and —2. 



ation of new branches balances the elimination of holes. 
The spectra exhibits two peaks; one at q « 1/R reflect- 
ing the typical distance of adjacent particles, and another 
characterizing the pattern size at = 1 j I* . In the quasi- 
stationary regime the characteristic length 1^,/R m 10 is 
independent of the system size, and only weakly depends 
on the particle density. 

The inscnsitivity of £* on particle density is in good 
agreement with the somewhat limited morphometric data 
available for the vasculature of quail embryos (see Fig. 3 
of [3j), where R « 20/im. LaRue et al characterized the 
vasculature with the average diameter of avascular ar- 
eas Ma and the average width of vascular segments My. 
From these two quantities, we estimate the characteristic 
pattern size as 1* « Ma + My and the volume fraction as 
cr « 1 — M 2 A j^. The five comparable data points cover 
volume fractions in the 0.2 < a < 0.9 range. The charac- 
teristic pattern size is = 80 ± 15 /^m with no obvious 
correlation between and a. 

As Fig. 4 demonstrates, the connectivity of the pat- 
tern depends on the particle density. We characterized 
the percolation transition by calculating the relative size 
of the largest interconnected cluster, P, and the volume 
fraction a by treating each particle as a disk of radius 
R and calculating the net space coverage of the configu- 
rations. The percolation threshold is at volume fraction 
a w 0.2: similar to the value reported in [lj|, and sub- 
stantially smaller than 0.67, the critical volume fraction 
for randomly placed overlapping disks [22| . 

To compare our model to that of [13], we analyze 
the structure of the critical cluster by calculating its 
mean density g as a function of radius r. Above the 
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FIG. 4: Network connectivity depends on particle density. 
Configurations of N = 2000 particles are shown in the sta- 
tionary phase for various values of N/L 2 particle density: 0.31 
(a), 0.44 (b), and 1.54 (c). The maximal cluster size (d) shows 
a percolation transition at a volume fraction below 0.2 (filled 
symbols), much less than the critical volume fraction of ran- 
domly placed disks (open symbols). The density-density au- 
tocorrelation of the critical cluster (e) is biphasic, and well 
approximated by g(r) ~ r -0,9 for r < r c and g(r) ~ r~ ' 2 for 
r > r c (solid lines). For comparison, the dashed line represent 
g(r) -r" - 5 . 

r/R ps 1 lower cut-off length g(r) is well approximated 
by a biphasic curve: g(r) ~ r -o.9±o.05 f or i < r < rc 
and g(r) ~ r -o.l7±o.i f or r ^ ^ r rpj^ crossover length 
r c ps QR is comparable with , the characteristic pattern 
size in the percolating regime. In agreement, for a wide 
range of particle densities, which includes a substantial 
regime above the percolation threshold, the S(q) at 
falls off to increasing q as S(q) ~ q^ 1 (Fig. 3a). While 
the obtained exponents are subject to finite size effects 
due to the limited scaling regime, both g(r) and S(q) is 
more compatible with linear structures below r c ps £ m) 
rather than the fractal-like behavior g(r) ~ r -0 5 seen in 
the chemoattractant model. 

In conclusion, we demonstrate that a preferential at- 
traction to elongated cells can be sufficient to explain the 
abundance of network-like structures in cell cultures and 
that it is likely to be an important component of vascu- 
logenic patterning. Network formation of mobile agents 
with spatially limited interaction range is a fundamental 



problem also occurring in technological fields: the es- 
tablishment of a self-organized communicating network 
between mobile robots is one recent example (23|. We 
expect that a similar sprouting mechanism in a system 
of self-propelled agents can create adaptive networks at 
low volume fractions. 
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